f(x) g(x)
x = 10^(-14): 0.4884981308350689 0.4999999999999988
x = 10^(-15): 0.44408920985006256 0.4999999999999999
x = 10^(-16): 0.0 0.529 Midterm Exam
Time Allowed: 1 hour 30 mins
Answer SECTION A and TWO of the three optional sections (B, C, D). If you have answered more than than the required two optional sections, you will be given credit for your best two optional sections. As a guide, I suggest you spend around 30 minutes on each section.
Please indicate your answers to multiple choice questions clearly with an ╳ to the left of your answers. Answer the optional sections on separate paper.
Please label any additional paper you use with your name and page number.
Electronic devices are not needed and not permitted in this examination. You may use your own notes, limited to 2 sides of letter sized paper. Answer the questions to Section A in the boxes provided. Additional paper may be requested.
29.1 A: [40 points] Compulsory Section
- [2 points] Suppose that a quantity x \not= 0 is approximated by \widetilde{x}. What is the relative error in this approximation?
\begin{align} \frac{|x - \widetilde{x}|}{|x|} \end{align}
- [3 points] Let x_n = \frac{1 + n^2}{ 2 n^3 } \cos n.
For which of the following sequences (\varepsilon_n) do we have x_n = O(\varepsilon_n) as n\to\infty?
Please select all that apply.
\varepsilon_n = n^{-2}
\varepsilon_n = n^{-1}
\varepsilon_n = 1
\varepsilon_n = n
Answer.
Since |x_n| \leq \frac{2 n^2}{ 2 n^3 } = \frac{1}{n}, we have x_n = O(n^{-1}) and thus also x_n = O(1) and x_n = O(n) as n \to \infty.
- [3 points] How many floating point numbers are there in the interval [1,2) with binary precision k = 3?
Answer.
Floating point numbers are of the form (-1)^s( 1 + f ) 2^{n} where f = 0 . b_1 b_2 \dots b_k in binary. Therefore, there are 2^k floating point numbers in the interval [1,2). When the binary precision is 3, there are 8 floating point numbers in the interval.
- [3 points] What is the relative condition number of the function f(x) = x - 1?
Answer.
\begin{align} \kappa_f(x) = \left| \frac{x f'(x)}{f(x)} \right| = \left| \frac{x}{x-1} \right| \end{align}
- [3 points] What scaling would you use in a graph to verify algebraic decay (i.e. O(n^{-\alpha}) for some \alpha > 0)?
Linear scale
Semi-logarithmic with the x-axis scaled
Semi-logarithmic with the y-axis scaled
Log-log plot (logarithmic scaling in both x- and y-axes)
Answer.
I would use a log-log scale. If y = n^{-\alpha}, then \log y = - \alpha \log n which is linear on a log-log plot.
- [3 points] Let us define
\begin{align*} f(x) = \frac{\sqrt{1+x}-1}{x}, \qquad g(x) = \frac{1}{\sqrt{1+x} + 1}. \end{align*}
One may verify that f and g are mathematically identical (but you don’t have to). The following table is the output of evaluating f and g on \{10^{-14}, 10^{-15}, 10^{-16}\} in double precision floating point. Why are the function values different?
Answer.
Subtractive cancellation. In the last line, \sqrt{1+x} = 1 in floating point arithmetic and so the result is 0.0.
- [4 points] Let f : [ 0, \tfrac\pi2 ] \to \mathbb R be given by f(x) = x^2\sin x - 1.
Show that there is a \xi \in [0,\tfrac\pi2] for which f(\xi) = 0.
Answer.
f is a continuous function with f(0) = -1 < 0 and f(\pi/2) = \pi^2/4 - 1 > 0. By the change of sign theorem, there exists some \xi \in [0,\pi/2] such that f(\xi) = 0.
f(π/2) = 1.4674011002723395
- [2 points] Here, we apply Newton’s method to f from Question 7.
Newton iteration on f(x) = x^2 sin(x) - 1 with x1 = π/2 ┌───────────┬──────────┬────────────────┬─────────────────────┐ │ iteration │ sequence │ residual │ approximate α │ │ n │ x[n] │ r[n]=|f(x[n])| │ log r[n]/log r[n-1] │ ├───────────┼──────────┼────────────────┼─────────────────────┤ │ 1.0 │ 1.5708 │ 1.4674 │ NaN │ │ 2.0 │ 1.10371 │ 0.0876848 │ -6.34694 │ │ 3.0 │ 1.06891 │ 0.00165227 │ 2.63171 │ │ 4.0 │ 1.06822 │ 6.52732e-7 │ 2.22338 │ │ 5.0 │ 1.06822 │ 1.01918e-13 │ 2.10043 │ └───────────┴──────────┴────────────────┴─────────────────────┘
Why do we use the residuals r_n := |f(x_n)| instead of the errors |x_n - \xi| when we are approximating the order of convergence (\alpha) in this method?
Answer.
\xi is unknown in advance.
Remark. Suppose that (x_n) \to \xi. Then, since f(x_n) = f(x_n) - f(\xi) = f'(\eta_n) (x_n - \xi) for some \eta_n between x_n and \xi, we have
\begin{align} \lim_{n\to\infty} \frac{ \log |f(x_{n+1})| }{ \log |f(x_n)| } % &= \lim_{n\to\infty}\frac{ \log |f'(\eta_{n+1})| + \log |x_{n+1} - \xi| }{ \log |f'(\eta_{n})| + \log |x_{n} - \xi| } \nonumber\\ % &= \lim_{n\to\infty}\frac{ \frac{\log |f'(\eta_{n+1})|}{\log |x_{n} - \xi|} + \frac{\log |x_{n+1} - \xi|}{\log |x_{n} - \xi|} }{ \frac{\log |f'(\eta_{n})|}{\log |x_{n} - \xi|} + 1 } \nonumber\\ % &= \lim_{n\to\infty} \frac{\log |x_{n+1} - \xi|}{\log |x_{n} - \xi|}. \end{align}
As a result, we may use the residuals to estimate the order of convergence.
- [4 points] Let g(x) := 1 + \frac{1}{x} and define x_1 = 1 and x_{n+1} = g(x_n). Suppose that (x_n) \to \xi.
What is the value of \xi?
Answer.
We have x_{n+1}\to \xi as n \to \infty. Since g is continuous, we have g(x_n)\to g(\xi) as n \to \infty. Therefore, \xi = 1+ \frac{1}\xi which is equivalent to \xi^2 - \xi - 1 = 0 and so \xi = \frac{1 \pm \sqrt{5}}2. If x_n > 0 then x_{n+1} >0. Since x_1 > 0 we have x_n > 0 for all n and thus \xi \geq 0. Therefore, \xi = \frac{1+\sqrt{5}}{2} (the golden ratio).
Remark.
It turns out the the sequence does indeed converge (can you find an interval on which g is a contraction and apply the contraction mapping theorem?)
a, b = 1.2, 2
g = x -> 1 + 1/x
plot( g, a, b, legend=false, lw =3 )
hline!( [a,b], linestyle=:dash)
x[n+1] = g(x[n]) with x_1 = 1.0 ┌───────────┬──────────┬────────────────┬─────────────────────┐ │ iteration │ sequence │ residual │ approximate α │ │ n │ x[n] │ r[n]=|f(x[n])| │ log r[n]/log r[n-1] │ ├───────────┼──────────┼────────────────┼─────────────────────┤ │ 1.0 │ 1.0 │ 1.0 │ NaN │ │ 2.0 │ 2.0 │ 0.5 │ -Inf │ │ 3.0 │ 1.5 │ 0.166667 │ 2.58496 │ │ 4.0 │ 1.66667 │ 0.0666667 │ 1.51139 │ │ 5.0 │ 1.6 │ 0.025 │ 1.36219 │ │ 6.0 │ 1.625 │ 0.00961538 │ 1.25902 │ │ 7.0 │ 1.61538 │ 0.003663 │ 1.20779 │ │ 8.0 │ 1.61905 │ 0.00140056 │ 1.17139 │ │ 9.0 │ 1.61765 │ 0.000534759 │ 1.14653 │ │ 10.0 │ 1.61818 │ 0.00020429 │ 1.12773 │ │ 11.0 │ 1.61798 │ 7.80275e-5 │ 1.11329 │ └───────────┴──────────┴────────────────┴─────────────────────┘ ξ = 1.618033988749895
Remark.
Notice that the first few terms of the sequence (x_n) are given by
\begin{align*} x_2 &= 1 + \frac{1}{x_1}, \qquad x_{3} = 1 + \frac{1}{x_2} = 1 + \frac{1}{1 + \frac{1}{x_1}}, \qquad x_{4} = 1 + \frac{1}{x_3} = 1 + \frac{1}{1 + \frac{1}{1 + \frac{1}{x_1}}}, \cdots. \end{align*}
As a result, you have given meaning to the infinite continued fraction
\begin{align*} \xi := 1 + \frac{1}{1 + \frac{1}{1 + \frac{1}{1 + \ddots}}} \end{align*}
- [4 points] Write down the Lagrange form of the polynomial of degree at most n interpolating f on X = \{ x_0, \dots, x_n \}.
Answer.
The Lagrange polynomials are \ell_j(x) := \prod_{0\leq k \leq n \, : \, k\not= j} \frac{x-x_k}{x_j - x_k} and the polynomial interpolation, p, of f on X is given by the following linear combination:
\begin{align} p(x) = \sum_{j=0}^n \ell_j(x) f(x_j) \end{align}
- [4 points] Consider the following finite difference table
\begin{align} x_0: \quad & f[x_0] = 0, \quad & f[x_0, x_1] = 1, \nonumber\\ x_1: \quad & f[x_1] = 0, \nonumber \end{align}
Which of the following statements are true?
Please select all that apply.
The table does not correspond to any function f,
The table determines a polynomial interpolant p of degreee at least 2,
We must have f'(x_1) = 1,
We must have x_0 = x_1,
Answer.
By the definition of divided differences, we have f(x_0) = f(x_1) = 0 and
\begin{align} 1 = f[x_0,x_1] = \begin{cases} \frac{f(x_0) - f(x_1)}{x_0 - x_1} &\text{if } x_0 \not= x_1 \\ f'(x_0) &\text{if } x_0 = x_1. \end{cases} \end{align}
- Therefore, we must have x_0 = x_1 and f'(x_0) = f'(x_1) = 1 (third and forth option should have been selected).
- This table only determines a polynomial, p, of degree at most 1 (and so the second option does not apply):
\begin{align} p(x) = f[x_0] + f[x_0,x_1] (x - x_0) = x - x_0. \end{align}
One can indeed show that p'(x_0) = 1 = f'(x_0) and p(x_0) = 0 = f(x_0).
- f(x) = x-x_0 is an example of a function with the divided difference table as shown (so the first option does not apply),
- [5 points] Write down the Trapezoid rule on [-1,1] and show that it is exact for all polynomials of degree at most 1.
Answer.
The trapezoid rule has weights w_0 = w_1 = 1 and nodes x_0 = -1 and x_1 = 1: i.e.
\begin{align} \int^{1}_{-1} f \approx f(-1) + f(1). \end{align}
When f(x) = \alpha, the integral is 2 \alpha = f(-1) + f(+1),
When f(x) = x, since f is odd, the integral is 0 = -1 + 1 = f(-1) + f(1).
Therefore, taking linear combinations, the rule is exact for all polynomials of degree 1.
Choose TWO sections from B, C, D (each worth 30 points)
29.2 B: Midpoint Rule
Suppose f is twice continuously differentiable. We will consider the following midpoint quadrature rule given by
\begin{align*} \int_a^b f(x) \approx (b-a) f\big( \tfrac{a+b}2 \big). \end{align*}
That is, there is only one node x_0 = \tfrac{a+b}{2} with corresponding weight w_0 = b-a.
- [3 points] Show that this midpoint rule is exact for all polynomials of degree at most 1,
Answer.
All polynomials of degree 1 may be written in the form f(x) = \alpha \big(x - \tfrac{a+b}{2}\big) + \beta. By anti-symmetry about \frac{a+b}2, we have
\begin{align} \int_a^b f &= \alpha \cdot 0 + \beta( b- a) \nonumber\\ % &= (b - a) \Big( \alpha \big(\tfrac{a+b}2 - \tfrac{a+b}{2}\big) + \beta \Big) \nonumber\\ % &= (b-a) f\big( \tfrac{a+b}{2} \big). \end{align}
- [3 points] Write down the polynomial p_1 of degree at most 1 with p_1\big( \tfrac{a+b}{2} \big) = f\big( \tfrac{a+b}{2} \big) and p_1'\big( \tfrac{a+b}{2} \big) = f'\big( \tfrac{a+b}{2} \big),
Answer.
\begin{align} p_1(x) = f\big( \tfrac{a+b}{2} \big) + f'\big( \tfrac{a+b}{2} \big) \big( x - \tfrac{a+b}2 \big) \end{align}
- [3 points] Show that there exists \xi_x \in [a,b] such that
\begin{align*} f(x) - p_1(x) = \frac{f''(\xi_x)}{2} \big(x - \tfrac{a+b}2\big)^2. \end{align*}
Answer.
By the Taylor remainder theorem, we have f(x) = f(\tfrac{a+b}2) + f'(\frac{a+b}2) (x - \xi) + \frac12 f''(\xi_x) (x - \tfrac{a+b}2)^2 = p_1(x) + \frac12 f''(\xi_x) (x - \tfrac{a+b}2)^2 for some \xi_x between x and \tfrac{a+b}2.
- [8 points] Show that
\begin{align*} \left| \int_{a}^b f(x) \mathrm{d}x - (b-a) f\big( \tfrac{a+b}{2} \big) \right| \leq \frac{\| f'' \|_{L^\infty([a,b])}}{24} (b-a)^3 . \end{align*}
Answer.
Since p_1\big( \tfrac{a+b}{2} \big) = f\big( \tfrac{a+b}{2} \big) and p_1 is a polynomial of degree at most 1 (and so the quadrature rule is exact for p_1), we have
\begin{align} \left| \int_{a}^b f(x) \mathrm{d}x - (b-a) f\big( \tfrac{a+b}{2} \big) \right| % &= \left| \int_{a}^b f(x) \mathrm{d}x - (b-a) p_1\big( \tfrac{a+b}{2} \big) \right| \nonumber\\ % &= \left| \int_{a}^b f(x) \mathrm{d}x - \int_a^b p_1(x) \mathrm{d}x \right| \nonumber\\ % &= \left| \int_{a}^b \left( f(x) - p_1(x) \right) \mathrm{d}x \right| \nonumber\\ % &= \left|\int_{a}^b \frac{f''(\xi_x)}{2} \big(x - \tfrac{a+b}2\big)^2 \mathrm{d}x\right| \nonumber\\ % &\leq \frac{\|f''\|_{L^\infty([a,b])}}{2} \int_{a}^b \big(x - \tfrac{a+b}2\big)^2 \mathrm{d}x \nonumber\\ % &= \frac{\|f''\|_{L^\infty([a,b])}}{6} \big(x - \tfrac{a+b}2\big)^3\Big|_a^b \nonumber\\ % &= \frac{\|f''\|_{L^\infty([a,b])}}{6} \frac{2(b-a)^3}{8}. \end{align}
Let z_{j} := a + j \frac{b-a}{n} for j = 0,\dots,n and note that \int_a^b f = \sum_{j=0}^{n-1} \int_{z_{j}}^{z_{j+1}} f.
- [5 points] Apply the midpoint rule on each interval [z_j, z_{j+1}] to write down a formula for the corrresponding composite rule on [a,b]. Give the weights \{w_j\} and nodes \{x_j\} of this rule.
Answer.
We find that
\begin{align} \int_a^b f &= \sum_{j=0}^{n-1} \int_{z_{j}}^{z_{j+1}} f \nonumber\\ % &\approx \sum_{j=0}^{n-1} \big( z_{j+1} - z_j \big) f\big( \tfrac{z_j + z_{j+1}}{2} \big) \nonumber\\ % &= \sum_{j=0}^{n-1} \frac{b-a}{n} f\big( a + (2j+1) \tfrac{b-a}{2n} \big) \end{align}
is a quadrature rule with nodes x_j = a + (2j+1) \frac{b-a}{2n} and weights w_j = \frac{b-a}n for j= 0,\dots,n-1.
- [4 points] Use the error estimate from Question 3 obtain an error estimate for the composite rule.
Answer.
Applying the error estimate on [z_j, z_{j+1}] gives an error of
\begin{align} \frac{\| f'' \|_{L^\infty([z_j,z_{j+1}])}}{24} (z_{j+1}-z_j)^3 % &\leq \frac{\| f'' \|_{L^\infty([a,b])}}{24} \left( \tfrac{b-a}{n} \right)^3 . \end{align}
Summing these errors for j = 0,\dots,n-1 gives
\begin{align} \left| \int_a^b f - \sum_{j=0}^{n-1} \frac{b-a}{n} f\big( a + (2j+1) \tfrac{b-a}{2n} \big) \right| % &\leq \frac{\| f'' \|_{L^\infty([a,b])}}{24} \frac{(b-a)^3}{n^2} . \end{align}
Remark. The rate of convergence of the composite rectangular rule was O(n^{-1}) and so, by choosing the midpoint instead of the left- (or right-) hand end-points, we get faster convergence at no extra computational cost.
- [4 points] Suppose that a = -1, b = 1. Show that the midpoint rule is exact when f is continuous and odd.
Answer.
Since f is odd, we have f(x) = -f(-x) and thus \int_{-1}^1 f = 0. Moreover, f is continuous, so f(0) = 0 and thus (1 - (-1)) f(\tfrac{-1 + 1}{2}) = 2 f(0) = 0. That is, the integral is equal to the quadrature rule in this case.
29.3 C: Simple Iteration
Define g(x) = 1 + \frac{1}{1 + x} and I := [1,2].
- [3 points] Show that g(x) \in I for all x \in I.
Answer.
Since g is decreasing, we have g(x) \in [ g(2), g(1) ] = [\tfrac{4}{3}, \tfrac{3}{2}] \subset [1,2] for all x \in I.
- [3 points] Which theorem guarantees the existence of a fixed point \xi = g(\xi) \in I.
Answer.
Since g is continuous, Brouwer’s fixed point theorem guarantees the existence of a fixed point \xi = g(\xi) \in I.
- [3 points] Compute the value of \xi.
Answer.
Since g is continuous, we have \xi = g(\xi) = 1 + \frac{1}{1 + \xi} which is equivalent to \xi^2 - 1 = 1 and thus \xi = \pm\sqrt{2}. Since x_n\in I for all n, we must have \xi \in I and so \xi = \sqrt{2}.
- [3 points] Show that L := \max_{x \in I} |g'(x)| = \frac14.
Answer.
We have g'(x) = - \frac{1}{(1+x)^2} which is also decreasing. As a result, we have |g'(x)| \leq \frac{1}{(1+1)^2} = \frac{1}{4} for all x \in I. This bound is attained by x = 1.
- [3 points] Explain why this means the iteration x_{n+1} = g(x_n) converges to \xi for all x_1 \in I.
Answer.
We have shown that g: I \to I is a contraction (|g(x) - g(y)| = |g'(z)| |x - y| \leq \frac{1}4|x-y| for some z \in I) and so by the contraction mapping theorem, we have x_{n+1} = g(x_n) converges to \xi for all x_1 \in I.
- [5 points] What is the order of convergence?
- [5 points] What is the asymptotic error constant?
Answer.
There exists \eta_n between x_n and \xi such that
\begin{align} \frac{|x_{n+1} - \xi|}{|x_n - \xi|} = \frac{|g(x_{n}) - g(\xi)|}{|x_n - \xi|} = | g'(\eta_n) |. \end{align}
As a result, taking the limit, we see that the order of convergence is 1 with asymptotic error constant
\begin{align} \mu = \lim_{n\to\infty}\frac{|x_{n+1} - \xi|}{|x_n - \xi|} = |g'(\xi)| = \frac{1}{(1+\sqrt{2})^2} \end{align}
- [5 points] Fix x_1 \in I. Find n for which the error bound |x_{n+1} - \xi| \leq \frac{1}{64} is guaranteed.
Answer.
We have |x_{n+1} - \xi| \leq \frac{1}4 |x_n - \xi| \leq \cdots \leq 4^{-n} |x_1 - \xi| \leq 4^{-n}. Since \frac{1}{64} = 4^{-3}, we would expect |x_{n+1} - \xi| \leq \frac{1}{64} whenever n \geq 3.
Remark.
Notice that the first few terms of the sequence (x_n) are given by
\begin{align*} x_2 &= 1 + \frac{1}{1 + x_1}, \qquad x_{3} = 1 + \frac{1}{1 + x_2} = 1 + \frac{1}{2 + \frac{1}{1 + x_1}}, \qquad x_{4} = 1 + \frac{1}{1 + x_3} = 1 + \frac{1}{2 + \frac{1}{2 + \frac{1}{1 + x_1}}}, \cdots. \end{align*}
As a result, you have given meaning to the infinite continued fraction
\begin{align*} \sqrt{2} = 1 + \frac{1}{2 + \frac{1}{2 + \frac{1}{2 + \ddots}}} \end{align*}
Remark.
While you didn’t need a computer to solve this problem, let’s just verify your findings…
x[n+1] = g(x[n]) with x_1 = 1.0 ┌───────────┬──────────┬───────────────┬─────────────────────┬───────────────┐ │ iteration │ sequence │ error │ approximate α │ μ (α = 1) │ │ n │ x[n] │ e[n]=|x[n]-ξ| │ log e[n]/log e[n-1] │ e[n]/e[n-1]^α │ ├───────────┼──────────┼───────────────┼─────────────────────┼───────────────┤ │ 1.0 │ 1.0 │ 0.414214 │ NaN │ NaN │ │ 2.0 │ 1.5 │ 0.0857864 │ 2.78644 │ 0.207107 │ │ 3.0 │ 1.4 │ 0.0142136 │ 1.73198 │ 0.165685 │ │ 4.0 │ 1.41667 │ 0.0024531 │ 1.41303 │ 0.172589 │ │ 5.0 │ 1.41379 │ 0.000420459 │ 1.29345 │ 0.171399 │ │ 6.0 │ 1.41429 │ 7.21519e-5 │ 1.22672 │ 0.171603 │ │ 7.0 │ 1.4142 │ 1.23789e-5 │ 1.18484 │ 0.171568 │ │ 8.0 │ 1.41422 │ 2.1239e-6 │ 1.156 │ 0.171574 │ │ 9.0 │ 1.41421 │ 3.64404e-7 │ 1.13495 │ 0.171573 │ │ 10.0 │ 1.41421 │ 6.25218e-8 │ 1.1189 │ 0.171573 │ │ 11.0 │ 1.41421 │ 1.0727e-8 │ 1.10627 │ 0.171573 │ │ 12.0 │ 1.41421 │ 1.84047e-9 │ 1.09606 │ 0.171573 │ └───────────┴──────────┴───────────────┴─────────────────────┴───────────────┘ theoretical value of μ = |g'(ξ)| = 0.1715728752538099 error after 3 iterations = 0.002453104293571595 < 1/64 = 0.015625
Choose TWO sections from B, C, D (each worth 30 points)
29.4 D: Newton’s Method for Repeated Roots
- [2 points] Show that f(x) = e^x - x - 1 satisfies f(0) = f'(0) = 0 and f''(0) \not= 0.
Answer.
We have f'(x) = e^x - 1 and f''(x) = f'''(x) = \cdots = e^x. Therefore, f(0) = e^0 - 0 - 1 = 0 and f'(0) = e^0 - 1 = 0 and f''(0) = e^0 = 1.
- [2 points] Here, we apply Newton’s method to f(x) = e^x - x - 1. What is the order of convergence?
┌───────────┬─────────────┬───────────────┬─────────────────────┐
│ iteration │ sequence │ error │ approximate α │
│ n │ x[n] │ e[n]=|x[n]-ξ| │ log e[n]/log e[n-1] │
├───────────┼─────────────┼───────────────┼─────────────────────┤
│ 1.0 │ 3.1 │ 3.1 │ NaN │
│ 2.0 │ 2.24624 │ 2.24624 │ 0.71527 │
│ 3.0 │ 1.512 │ 1.512 │ 0.51088 │
│ 4.0 │ 0.939627 │ 0.939627 │ -0.150621 │
│ 5.0 │ 0.542328 │ 0.542328 │ 9.826 │
│ 6.0 │ 0.295555 │ 0.295555 │ 1.99205 │
│ 7.0 │ 0.155046 │ 0.155046 │ 1.52927 │
│ 8.0 │ 0.0795256 │ 0.0795256 │ 1.35817 │
... ... ... ...
│ 12.0 │ 0.0050958 │ 0.0050958 │ 1.15071 │
│ 13.0 │ 0.00255006 │ 0.00255006 │ 1.13113 │
│ 14.0 │ 0.00127557 │ 0.00127557 │ 1.116 │
│ 15.0 │ 0.000637923 │ 0.000637923 │ 1.10398 │
│ 16.0 │ 0.000318995 │ 0.000318995 │ 1.0942 │
│ 17.0 │ 0.000159506 │ 0.000159506 │ 1.0861 │
│ 18.0 │ 7.97552e-5 │ 7.97552e-5 │ 1.07927 │
│ 19.0 │ 3.98781e-5 │ 3.98781e-5 │ 1.07345 │
│ 20.0 │ 1.99392e-5 │ 1.99392e-5 │ 1.06843 │
│ 21.0 │ 9.96963e-6 │ 9.96963e-6 │ 1.06404 │
└───────────┴─────────────┴───────────────┴─────────────────────┘
Answer.
The computational order of convergence appears to approach 1.
Consider the modified Newton’s method
\begin{align*} x_{n+1} = x_n - 2 \frac{f(x_n)}{f'(x_n)}. \end{align*}
- [2 points] We apply this method to f(x) = e^x - x - 1. Compare the table to that of the standard Newton method.
┌───────────┬─────────────┬───────────────┬─────────────────────┐ │ iteration │ sequence │ error │ approximate α │ │ n │ x[n] │ e[n]=|x[n]-ξ| │ log e[n]/log e[n-1] │ ├───────────┼─────────────┼───────────────┼─────────────────────┤ │ 1.0 │ 3.1 │ 3.1 │ NaN │ │ 2.0 │ 1.39248 │ 1.39248 │ 0.292634 │ │ 3.0 │ 0.313183 │ 0.313183 │ -3.50653 │ │ 4.0 │ 0.0163206 │ 0.0163206 │ 3.54474 │ │ 5.0 │ 4.43937e-5 │ 4.43937e-5 │ 2.43539 │ │ 6.0 │ 3.27151e-10 │ 3.27151e-10 │ 2.17918 │ └───────────┴─────────────┴───────────────┴─────────────────────┘
Answer.
We have convergence in fewer iterations (21 vs 6) and the computational order of convergence seems to be around 2.
Now suppose that f:[a,b]\to\mathbb R is four times continuously differentiable with root \xi \in (a,b) of multiplicity 2: that is, f(\xi) = f'(\xi) = 0 and f''(\xi) \not= 0. Further, we suppose that there exists A > 0 such that
\begin{align*} \left| \frac{f'''(x) + f^{(4)}(x) (x - \xi)}{f''(y)} \right| \leq A \end{align*}
for all x,y \in [a,b]. The aim of this section is to show that the modified Newton method converges at least quadratically to \xi for all x_1 \in [a,b] with |x_1 - \xi| \leq A^{-1}.
Remark.
Recall that in lectures, we showed a similar result but for the standard Newton method under the condition that \left| \frac{f''(x)}{f'(y)} \right| remains bounded for x and y near \xi.
Notice that the modified Newton method can be written
\begin{align*} x_{n+1} - \xi &= \frac{ f'(x_n)(x_n - \xi) - 2 f(x_n) }{ f'(x_n) }. \end{align*}
- [8 points] By expanding the numerator g(x) := f'(x) (x - \xi) - 2 f(x) and denominator h(x) := f'(x), show that there exist \eta_n, \theta_n between x_n and \xi such that
\begin{align*} x_{n+1} - \xi &= \frac {\frac1{3!} g'''(\eta_n) (x_n - \xi)^3} {h'(\theta_n) (x_n - \xi) } % = \frac16 \left[ \frac {f'''(\eta_n) + f^{(4)}(\eta_n) (\eta_n - \xi) } {f''(\theta_n)} \right] (x_n - \xi)^2. \end{align*}
Answer.
We take derivatives of g:
\begin{align} g(x) &= f'(x) (x - \xi) - 2 f(x)\nonumber\\ g'(x) &= f''(x) (x - \xi) - f'(x) \nonumber\\ g''(x) &= f'''(x) (x - \xi) \nonumber\\ g'''(x) &= f^{(4)}(x) (x - \xi) + f'''(x). \end{align}
Notice that g(\xi) = g'(\xi) = g''(\xi) = 0 and so there exists \eta_n between x_n and \xi such that
\begin{align} g(x_n) &= g(\xi) + g'(\xi) (x_n - \xi) + \frac{1}2 g''(\xi) (x_n - \xi)^2 + \frac{1}{6} g'''(\eta_n) (x_n - \xi)^3 \nonumber\\ % &= \frac{1}{6} \Big[ f'''(\eta_n) + f^{(4)}(\eta_n) (\eta_n - \xi) \Big] (x_n-\xi)^3. \end{align}
Moreover, there exists \theta_n such that f'(x_n) = f'(\xi) + f''(\theta_n) (x_n - \xi) = f''(\theta_n) (x_n - \xi) and so
\begin{align} x_{n+1} - \xi &= \frac16 \left[ \frac {f'''(\eta_n) + f^{(4)}(\eta_n) (\eta_n - \xi) } {f''(\theta_n)} \right] (x_n - \xi)^2 \end{align}
- [4 points] Hence show that, if |x_1 - \xi| \leq A^{-1}, then |x_{n+1} - \xi| \leq 6^{-n} |x_1 - \xi| and thus the modified Newton’s method converges.
Answer.
By assumption,
\begin{align} \left| \frac {f'''(\eta_n) + f^{(4)}(\eta_n) (\eta_n - \xi)} {f''(\theta_n)} \right| \leq A, \end{align}
and so, if |x_n - \xi| \leq A^{-1}, then
\begin{align} |x_{n+1} - \xi| \leq \frac16 A |x_n - \xi| |x_n - \xi| % \leq \frac16 |x_{n} - \xi|. \end{align}
In particular, |x_{n+1} - \xi| \leq A^{-1} and we may apply the bound inductively to obtain |x_{n+1} - \xi| \leq 6^{-n} |x_1 - \xi| and so this modified Newton’s method converges.
- [4 points] Show that (x_n) converges at least quadratically to \xi.
Answer.
Since \eta_n, \theta_n are between x_n and \xi, and (x_n) \to \xi, we have
\begin{align} \lim_{n\to\infty} \frac{|x_{n+1} - \xi|}{|x_n - \xi|^2} &= \frac16 \left| \frac {f'''(\eta_n) + f^{(4)}(\eta_n) (\eta_n - \xi) } {f''(\theta_n)} \right| \nonumber\\ % &= \frac16 \left| \frac {f'''(\xi)} {f''(\xi)} \right| \end{align}
- [4 points] What is the theoretical value of the asymptotic error constant when the modified Newton method is applied to f(x) = e^x - x - 1?
Answer.
Recall that, in this case, \xi = 0, and f'(x) = e^x - 1, f''(x) = f'''(x) = e^x. Therefore, the asymptotic error constant is \frac16 \frac{e^0}{e^0} = \frac16.
- [4 points] Show that the modified Newton method converges in 1 iteration for functions of the form f(x) = (x - \xi)^2.
Answer.
Note that f(x)= (x - \xi)^2 has a root of multiplicity 2 at \xi. We have x_{n+1} = x_n - 2 \frac{ (x_n - \xi)^2 }{ 2(x_n - \xi)} = x_n - (x_n - \xi) = \xi so that x_{n+1} = \xi for all n \geq 1 regardless of x_1 \in \mathbb R.
Remark.
In general, if \xi is a root of f of multiplicity m, one can define the following modified Newton method:
\begin{align} x_{n+1} = x_n - m \frac{f(x_n)}{f'(x_n)}. \end{align}
Then, if f is sufficiently smooth and x_1 is sufficiently close to \xi, the modified Newton method converges quadratically with asymptotic error constant
\begin{align} \mu = \frac{1}{m(m+1)} \left|\frac{f^{(m+1)}(\xi)}{f^{(m)}(\xi)} \right| \end{align}
We quickly verify this numerically for the function f(x) = e^x - \frac12 x^2 - x - 1 which has root \xi = 0 of multiplicity 3.
x[n+1] = g(x[n]) with x_1 = 1.0 ┌───────────┬──────────┬───────────────┬─────────────────────┬───────────────┐ │ iteration │ sequence │ error │ approximate α │ μ (α = 1) │ │ n │ x[n] │ e[n]=|x[n]-ξ| │ log e[n]/log e[n-1] │ e[n]/e[n-1]^α │ ├───────────┼──────────┼───────────────┼─────────────────────┼───────────────┤ │ 1.0 │ 1.0 │ 0.414214 │ NaN │ NaN │ │ 2.0 │ 1.5 │ 0.0857864 │ 2.78644 │ 0.207107 │ │ 3.0 │ 1.4 │ 0.0142136 │ 1.73198 │ 0.165685 │ │ 4.0 │ 1.41667 │ 0.0024531 │ 1.41303 │ 0.172589 │ │ 5.0 │ 1.41379 │ 0.000420459 │ 1.29345 │ 0.171399 │ │ 6.0 │ 1.41429 │ 7.21519e-5 │ 1.22672 │ 0.171603 │ │ 7.0 │ 1.4142 │ 1.23789e-5 │ 1.18484 │ 0.171568 │ │ 8.0 │ 1.41422 │ 2.1239e-6 │ 1.156 │ 0.171574 │ │ 9.0 │ 1.41421 │ 3.64404e-7 │ 1.13495 │ 0.171573 │ │ 10.0 │ 1.41421 │ 6.25218e-8 │ 1.1189 │ 0.171573 │ │ 11.0 │ 1.41421 │ 1.0727e-8 │ 1.10627 │ 0.171573 │ │ 12.0 │ 1.41421 │ 1.84047e-9 │ 1.09606 │ 0.171573 │ └───────────┴──────────┴───────────────┴─────────────────────┴───────────────┘ theoretical value of μ = |g'(ξ)| = 0.1715728752538099 error after 3 iterations = 0.002453104293571595 < 1/64 = 0.015625